Role of Sr doping and external strain on relieving bottleneck of oxygen diffusion in La2−xSrxCuO4−δ

In many complex oxides, the oxygen vacancy formation is a promising route to modify the material properties such as a superconductivity and an oxygen diffusivity. Cation substitutions and external strain have been utilized to control the concentration and diffusion of oxygen vacancies, but the mechanisms behind the controls are not fully understood. Using first-principles calculations, we find how Sr doping and external strain greatly enhances the diffusivity of oxygen vacancies in La2−xSrxCuO4−δ (LSCO) in the atomic level. In hole-doped case (2x > δ), the formation energy of an apical vacancy in the LaO layer is larger than its equatorial counterpart by 0.2 eV that the bottleneck of diffusion process is for oxygen vacancies to escape equatorial sites. Such an energy difference can be reduced and even reversed by either small strain (< 1.5%) or short-range attraction between Sr and oxygen vacancy, and in turn, the oxygen diffusivity is greatly enhanced. For fully compensated hole case (2x ≦ δ), the formation energy of an apical vacancy becomes too high that most oxygen vacancies cannot move but would be trapped at equatorial sites. From our electronic structure analysis, we found that the contrasting change in the formation energy by Sr doping and external strain is originated from the different localization natures of electron carrier from both types of oxygen vacancies.

Cation substitutions and external strain has been utilized to control the electronic and diffusion properties of oxygen vacancies, but the mechanisms behind the controls are not fully understood. A part of the reason is that though the effects of oxygen vacancy strongly depend on crystallographic location (apical or equatorial) in its bulk structures, the measurements are performed in indirect ways such as Raman scattering 25 , lattice expansion 20 , and oxygen tracer diffusion 29 . In this work, to directly determine the crystallographic locations of oxygen vacancies in La 2−x Sr x CuO 4 (LSCO) and the role of them in the electronic and diffusion properties under external strain and doping, we scrutinize the energetic and electronic properties of oxygen vacancies using first-principles calculations. We found that the formation energies ( E V ) of apical ( aV O ) and equatorial oxygen vacancies ( eV O ) depend strongly on doping concentrations and strain, but in different manners. Implications on the oxygen diffusivity and how external strain can effectively remove the bottleneck of oxygen diffusion in bulk LSCO will be discussed. From the vacancy level analysis, the contrasting change in the formation energy by Sr doping and external strain are from the different localization natures of electron carrier from apical and equatorial vacancy. This also indicate that it is nearly impossible to achieve electron doping by introducing oxygen vacancies.

Theoretical approaches
We performed first-principles calculations using density functional theory (DFT) 30 as implemented in Vienna ab initio simulation package (VASP) 31,32 . The spin-polarized generalized gradient approximation (SGGA) with Perdew-Burke-Ernzerhof (PBE) functional 33 was used to treat the exchange and correlation functional. The electron wave functions were expanded by plane-wave basis with a kinetic energy cutoff of 400 eV. To introduce a low-temperature structure with alternatingly tilted octahedra (See "Atomic structure of La 2 CuO 4 " section and Table 1 for a detailed description of the structure), we used a 2 √ 2 × 2 √ 2 R45 • × 2 supercell containing 16 formula units and a 4 × 4 × 1 k-point mesh that uniformly samples the Brillouin zone (BZ) of the supercell. To obtain accurate densities of states, the BZ were sampled with a denser 8 × 8 × 4 k-point grid. All supercell structures were optimized until the force exerting on every atom became less than 0.02 eV/Å. The electron correlation effect associated with the Cu 3d orbitals was incorporated by the on-site Coulomb repulsion 34 [35][36][37] . We choose U eff = 7.0 eV, which reproduces experimental Cu magnetic moment and band gap (See Table 1). Although E V itself depends on U eff , the relative E V between aV O and eV O is much less sensitive on the choice of U eff . For unstrained LSCO, we checked this using U eff = 4.0 eV, which is from the reaction enthalpy of copper oxides 38 , confirming that our results on competitive formation of aV O and eV O does not depends on the choice of U eff (see Supplementary Information Section 1). For the convenience of analysis, the lattice constant along the c-axis was fixed while applying the biaxial strain along the a-and b-axes. The lattice constant along the c-axis is affected by elastic moduli, Sr doping 39 and oxygen vacancy concentration 20 , and the latter two effects are negligible in the vacancy formation energy. From our calculations, under the biaxial strain of ε (|ε|≤ 1.5%), c-axis is strained by -0.7 ε. Ignoring this elastic modulus effect slightly changes E V in strained LSCO, but this amount is quite small (by definition, zero for unstrained LSCO) that our results do not be affected by the fixed c-axis calculations. (See Supplementary Information Section 1).

Results and discussion
Atomic structure of La 2 CuO 4 . La 2 CuO 4 can be understood as a stacked structure of CuO 2 layer (A) and LaO layers (B), with (AB)B(AB)B··· layers repeating as shown in Fig. 1a. As the temperature increases, La 2 CuO 4 undergoes a successive phase transition from the low-temperature octahedral (LTO) phase to the low-temperature tetragonal (LTT) phase and then to the high-temperature tetragonal (HTT) phase. In the LTO and LTT phases, CuO 6 octahedra are alternately tilted along the [100] and [110] axes (Fig. 1b), respectively, while they are presumably randomly tilted in the HTT phase. The total energy of the LTO phase is lower than the LTT phase by 4.0-4.4 meV per formula unit 40,43 , and the small energy difference can be easily reversed by cation doping 41 . Assuming the choice of a specific low-energy phase has little effect on the formation energy of oxygen vacancies, we adopt the LTT phase in our calculations. We constructed an optimized supercell structure of the LTT phase consisting of 16 formula units, as shown in Fig. 1a, with alternatingly-tilted octahedra (see Fig. 1c). In the LTT phase, there are three inequivalent oxygen sites: one apical site ( aV O ) and two equatorial sites ( eV O , and eV * O ) , with no and an out-of-plane displacement, respectively (see Fig. 1b). Table 1   www.nature.com/scientificreports/ energy difference is 28 meV, which is quite weak for Sr atoms to form any type of ordering during the hightemperature process of synthesis. On the other hands, between Sr atom and oxygen vacancy, much larger attractive interaction exists. At zero strain with n = 2, E V of aV O is 1.98 eV when Sr is far away for aV O (6.09 Å), and it becomes lower to 1.76 eV when Sr is at the nearest location with aV O (2.59 Å). For eV O , the attraction is smaller and short-ranged that E V changes from 1.77 eV to 1.71 eV only when Sr atom is the nearest neighbor of eV O (2.70 Å). More details can be found in Supplementary Information Section 3. The strength of Sr-vacancy interaction is strong enough to affect the competitive formation of eV O and aV O at hole-doped sample as we will see soon.  (Fig. 2b,c) indicate the effect of octahedral tilt is negligible and we will focus only on eV * O from now on. Sr atoms provide n holes per supercell and they recombine with the oxygen vacancy making it a charged vacancy. As n increases, regardless of oxygen sites and strain, E V greatly decreases until n = 2, while further Sr doping beyond n = 2 scarcely change E V , indicating the stability of + 2 charged vacancy ( V +2 O ). Interestingly, the energy difference between eV * O and aV O in 0% strain case significantly reduces from 1.84 eV in the undoped sample to 0.20 eV in the sample with the + 2 charge states. The general trend that increasing Sr doping greatly reduces E V until n = 2 can be understood as follows. Once an oxygen vacancy forms, formally two electrons are released and they are localized around the vacancy. They now lose the energy gain from the large electron affinity of oxygen that E V is strongly affected by where and how tightly these two electrons are localized. In this sense, a Sr dopant (substituting La) can be thought to provides a location with a large electron affinity for one electron. The reason why the reduction of E V on Sr doping is so www.nature.com/scientificreports/ different for aV O and eV O can be explained with the different localization natures of electrons, and the detailed analysis will be presented in "Charge distribution and structural distortion by oxygen vacancies" section.
Our results imply that most oxygen vacancies in the undoped sample are located at equatorial sites. In the view of diffusion process via the vacancy mechanism, the bottleneck for the bulk diffusion is the equatorial-to-apical sites hopping, which is the only channel for diffusion along the c-axis. This hopping is nearly impossible to occur at the operation temperature of practical solid oxide fuel cells (< 1000 °C) because the lower bound of activation barrier is the relative E V between aV O and eV O (> 1 eV). However, in the hole-doped sample, oxygen vacancy can hop between equatorial and apical sites due to the reduced formation energy difference, and for zero strain and n = 2 case, calculated energy barrier between apical-to-equatorial site is 0.52 eV. The attractive interaction between Sr and aV O reduces E V of aV O , and in turn, it further reduces the energy barrier up to ~ 0.1 eV when a Sr atom is at the nearest neighbor of aV O . Experimentally, the formation of aV O in hole-doped samples was evidenced by confocal Raman microscopy 25 . Though the signal is interpreted as selective formation of apical vacancies, we think further investigations are called for to confirm it because from our calculation, the energy difference between both types of vacancies is not so large for the prevalence of a specific type.
The effect of strain on E V is quite different for eV * O and aV O . In the case of aV O , E V decreases linearly with a tensile strain at a similar rate for different Sr-doped samples. On the contrary,E V of eV * O shows relatively weak dependence on the strain and even non-monotonic behavior for neutral and + 1 charged vacancies ( V +1 O ). This difference is mainly from the distinctive bonding nature of apical and equatorial oxygen atom where the former is more ionic. We will investigate this further in "Charge distribution and structural distortion by oxygen vacancies" section by analyzing the structural distortions and accompanying internal stress from the oxygen vacancies. The small energy difference between + 2 charged aV O and eV * O is further reduced when biaxial strain is applied and even reversed at about 1-1.5% of strain as shown in Fig. 2d, which is achievable for epitaxially grown samples. For La 1.85 Sr 0.15 CuO 4 , it was clearly demonstrated that the tensile strain greatly accelerates reversible redox reaction 20 . Because the amount of reduction in E V on the biaxial strain is not so drastic for explaining the experiments, the acceleration was attributed to the surface exchange kinetics 20 . But for the whole sample to be uniformly doped, the rate-determining step should be the bulk diffusion including equatorial-to-apical hopping process. Our calculation indicates that the major role of tensile strains is to reduce E V of aV O , which does in turn, make the hopping process easier.

Charge distribution and structural distortion by oxygen vacancies.
Once an oxygen vacancy forms, formally two electrons are released and they are weakly bound at (or transferred to) nearby Cu atoms having d 9 configurations. This can be quantified using Bader charge analysis 48 where total electron charge is divided and assigned to each atom according to Bader's zero flux surface scheme. The assigned charge is called Bader charge and useful for tracing the charge transfer of a system. For eV * O , 1.15 electrons are transferred from oxygen atoms to Cu atoms, and 86% of them uniformly occupy the vacant d-orbitals of the two nearest Cu atoms, as shown in Fig. 3a. This results in a complete suppression of the magnetic moments of those two Cu atoms. Figure 3b shows charge allocations to the two Cu atoms for different charge states. For + 1 charge state, due to the strong on-site Coulomb repulsion of the Cu d-orbitals, only one of two copper atoms loses its electron, breaking the symmetric charge distribution. For + 2 charge states, two Cu atoms return to original d 9 configurations, recovering their magnetic moment (0.579 µ B ) almost to the pristine bulk value (0.603 µ B ). The structural distortion caused by eV * O is mainly from a weak attraction between the two Cu atoms (Cu1 and Cu2). As shown in Fig. 3c, the weak attraction causes two Cu atoms to approach each other and simultaneously rotate www.nature.com/scientificreports/ the neighboring octahedra, indicated by the blue arrows. The interatomic distance between the two Cu atoms is smallest for neutral eV * O (3.511 Å) and gradually recovers to the pristine bulk value (3.832 Å) as being positively charged (specifically, 3.757 Å for eV +2 O ). For aV O , one of the two released electrons from the oxygen vacancy is captured by Cu3 atom at the bottom of the vacancy and the remaining one is transferred to the lower CuO 2 plane and spatially extended as shown in Fig. 3d. This is somewhat contrasted to F-center in an ionic crystal where the released electrons from an anionic vacancy are well localized at the vacancy. Within Bader charge analysis (Fig. 3e), 1.51 electrons are transferred from the oxygen atom and only a third of them occupy the hole state of the bottom Cu3 atom in no charge state. For + 1 charge state, the extended states first recombine with the hole, whereas for + 2 charge states, Cu3 recovers its magnetic moment, 0.580 µ B . This charge redistribution causes copper layers and octahedrons to be gradually rearranged and makes the structure unstable. The absence of the negative oxygen ion causes the four positive lanthanum ions to repel each other by 0.334 Å, as shown in Fig. 3f. Similar with eV O , vacancy charging mitigates the structural distortions around oxygen vacancies.
As shown in Fig. 2a, E V of aV O decrease linearly with tensile strain at a similar rate regardless of their charge states. In most ionic crystals and perovskite materials such as BaTiO 3 , SrTiO 3 , SrZrO 3 , and PbZrO 3 , cations around anionic vacancy move away from each other 49 , indicating the strong Coulomb repulsions between the cations. This induces strong local compressive stress in the crystal and if the crystal is allowed to relax, the equilibrium lattice constants along the a-and b-axes increase by 0.6% for n = 0 case of our supercell geometry. The dependence on charge states is minute since the electron carriers from aV O are away from the apical sites that they play little role in reducing the Coloumb repulsion energy between the La atoms. Becase the tensile strain can mitigate the local compressive stress, this explains the decrease in E V of aV O on tensile strain at a similar rate regardless of their charge states.
On the other hand, E V 's of eV * O (or eV O ) show a relatively weak, non-monotonic dependence on both strain and charge states, as shown in Fig. 2c. In CuO 2 plane, the attraction between two Cu atoms around eV * O induces local tensile stress while in the neighboring two LaO planes, small but finite compressive stress is present, which we can infer from the displacement of La atoms. The non-monotonic behavior comes from the competition of the two opposite stresses. As eV * O is positively charging, the attraction becomes weaker and for + 2 charge states, compressive stress in LaO planes prevails and this explains weak linear decrease of E V on tensile biaxial strain. Figure 4a shows the total density of states (TDOS) and projected density of states (PDOS) of La 2 CuO 4 . From orbital-projected DOS analysis (not shown here), it was found that the fully occupied oxygen 2p orbitals are hybridized with the Cu 3d orbitals located near the top of valence band, while the d x 2 −y 2 hole states of Cu form the bottom of the conduction band. In Fig. 4b, eV * O induces two vacancy states (Cu1 and Cu2 states) at 0.8 eV above the valence band maximum. The energy level of the d 9 states of the two nearest Cu atoms from the vacancy is also changed due to the breaking of the Cu-O bond, which are manifested as the four in-gap states indicated by a red rectangle in Fig. 4b. The distribution of two vacancy states are plotted in Fig. 4d. It clearly shows that the additional electrons are localized in the Cu atoms. If eV * O is singly charged, one of the Cu states becomes unoccupied and the energy degeneracy of d 9 states are slightly lifted as shown in Fig. 4b. As explained in "Formation energy of oxygen vacancies as a function of Sr doping and biaxial strain", if the sample is hole-doped and oxygen vacancy is less than half of the hole concen- www.nature.com/scientificreports/ tration, most electrons in the vacancy states recombine with two holes leaving two unoccupied in-gap states, as in n = 2 case. In the case of aV O , one of two released electron forms a vacancy state associated with Cu3 at 1.3 eV above the valence band maximum while the other electron does a shallow level which is very close to the conduction band minimum as shown in Fig. 4c. Though aV O seems to play an electron donor, E V of it is very high compared to that of eV * O indicating that aV O is hard to form, or fast diffuses into the equatorial sites. In other words, oxygen vacancy in LSCO can compensate the holes in the sample but electron doping by overcompensation ( 2δ > x ) may not occur. Charge distribution of the vacancy state (Cu3) is shown in Fig. 4d. Other in-gap states in down spin channel come from one of the Cu3 d 9 states due to Cu-O bond breaking (blue dotted arrow). If aV O is singly charged, the electron in the shallow level first recombines with the hole. In the + 2 charge state, unlike eV * O , hole recombination makes Cu3 state to merge into conduction band, leaving no in-gap states.

Conclusion
We have calculated E V of eV * O and aV O in LSCO considering their charge states and strain dependence. In holedoped cases, most oxygen vacancies become + 2 charge states and compensate the hole density regardless of crystallographic location. E V of eV * O is lower than that of aV O by 0.2 eV, but the small energy difference can be easily overcome by 1.5% of tensile strain or Sr-aV O attractions, enhancing the oxygen diffusivity in Sr-overdoped samples. Even though the number of vacancies is sufficiently large to fully compensate the holes, the sample is hard to be electron-doped because the electrons from the equatorial vacancies are strongly bound around the vacancies. Neutral aV O can be an electron donor due to shallow vacancy level, but E V of it is too high to form